Split into three parts:
Projection parameters
Loads data frame with population parameters and explores the structure
of stage-based matrix models.
Projection models
Run deterministic and stochastic projections, combine results and export
for future use.
Model lookup
Summarise projection models and export for future use in scenario
modelling.
library(plyr)
library(tidyverse)
library(popdemo)
library(popbio)
library(Rcompadre)
library(Rage)
Load functions.
source('R/pop01_param_poun.R')
source('R/pop03_doproj.R')
source('R/pop03_doproj_stoch.R')
Each row holds a unique set of parameters for a matrix population model. This reflects diverse recovery (conservation action) and extinction (threat) scenarios.
dt <- pop01_param_poun()
nf <- 10
dt$adultF_n <- nf
ceiling_threshold <- nf + (nf * 0.2)
The first few rows should look like this….
| species | type | first_year | a1 | a2 | a3 | a4 | b1 | b2 | b3 | b4 | c1 | c2 | c3 | c4 | d1 | d2 | d3 | d4 | akey | adultF_n |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Podocnemis unifilis | base | 0.0 | 0 | 0 | 0 | 10.74 | 0.0 | 0.33 | 0 | 0 | 0 | 0.17 | 0.29 | 0 | 0 | 0 | 0.11 | 0.93 | poun_base_0 | NA |
| Podocnemis unifilis | base | 0.1 | 0 | 0 | 0 | 10.74 | 0.1 | 0.33 | 0 | 0 | 0 | 0.17 | 0.29 | 0 | 0 | 0 | 0.11 | 0.93 | poun_base_0.1 | NA |
| Podocnemis unifilis | base | 0.2 | 0 | 0 | 0 | 10.74 | 0.2 | 0.33 | 0 | 0 | 0 | 0.17 | 0.29 | 0 | 0 | 0 | 0.11 | 0.93 | poun_base_0.2 | NA |
Each row can be easily converted for use in a stage-based matrix population model. The following code shows how the first row becomes a stage-based matrix.
stage_names <- c("a1", "a2", "a3", "a4",
"b1", "b2", "b3", "b4",
"c1", "c2", "c3", "c4",
"d1", "d2", "d3", "d4")
vpop <- unlist(dt[3 , stage_names])
pop_mat <- matrix(vpop, byrow = TRUE, ncol=4)
dimnames(pop_mat) <- list(c("a", "b", "c", "d"),
c( "a", "b", "c", "d"))
The matrix should look like this:
| a | b | c | d | |
|---|---|---|---|---|
| a | 0.000000 | 0.000000 | 0.000000 | 10.741500 |
| b | 0.200000 | 0.333333 | 0.000000 | 0.000000 |
| c | 0.000000 | 0.166667 | 0.285714 | 0.000000 |
| d | 0.000000 | 0.000000 | 0.114286 | 0.930000 |
This matrix includes information of population growth and survival and reproduction.
See https://cran.r-project.org/web/packages/Rage/vignettes/a01_GettingStarted.html
We can show the different matrix components converting to the structure recommended by …
meta <- data.frame(idNum = 3,
SpeciesAccepted = dt[3 , 'species'],
type = dt[3 , 'type'],
first_year = dt[3 , 'first_year'])
stageInfo <- list(
data.frame(
MatrixClassOrganized = rep("active", 4),
MatrixClassAuthor = c("eggs/hatchling", "juvenile-early", "juvenile-late", "adult")
))
# Simple full A matrix
l_pop_mat <- list(mpm = pop_mat)
x <- cdb_build_cdb(mat_a = l_pop_mat, metadata = meta, stages = stageInfo)
# Seperate into U, F, C
mat_u1 <- rbind(
c(dt[3, 'a1'], dt[3, 'a2'], dt[3, 'a3'], 0),
c(dt[3, 'b1'], dt[3, 'b2'], dt[3, 'b3'], dt[3, 'b4']),
c(dt[3, 'c1'], dt[3, 'c2'], dt[3, 'c3'], dt[3, 'c4']),
c(dt[3, 'd1'], dt[3, 'd2'], dt[3, 'd3'], dt[3, 'd4'])
)
mat_f1 <- rbind(
c(0.0, 0.0, 0.0, dt[3, 'a4']),
c(0.0, 0.0, 0.0, 0.0),
c(0.0, 0.0, 0.0, 0.0),
c(0.0, 0.0, 0.0, 0.0)
)
mat_c1 <- rbind(
c(0.0, 0.0, 0.0, 0.0),
c(0.0, 0.0, 0.0, 0.0),
c(0.0, 0.0, 0.0, 0.0),
c(0.0, 0.0, 0.0, 0.0)
)
l_u <- list(m_u = mat_u1)
l_f <- list(m_f = mat_f1)
l_c <- list(m_c = mat_c1)
my_comadre <- cdb_build_cdb(
mat_u = l_u, mat_f = l_f, mat_c = l_c,
metadata = meta,
stages = stageInfo
)
matA(my_comadre)
## [[1]]
## [,1] [,2] [,3] [,4]
## [1,] 0.0 0.000000 0.000000 10.7415
## [2,] 0.2 0.333333 0.000000 0.0000
## [3,] 0.0 0.166667 0.285714 0.0000
## [4,] 0.0 0.000000 0.114286 0.9300
matA(x)
## [[1]]
## a b c d
## a 0.0 0.000000 0.000000 10.7415
## b 0.2 0.333333 0.000000 0.0000
## c 0.0 0.166667 0.285714 0.0000
## d 0.0 0.000000 0.114286 0.9300
matU(my_comadre)
## [[1]]
## [,1] [,2] [,3] [,4]
## [1,] 0.0 0.000000 0.000000 0.00
## [2,] 0.2 0.333333 0.000000 0.00
## [3,] 0.0 0.166667 0.285714 0.00
## [4,] 0.0 0.000000 0.114286 0.93
matF(my_comadre)
## [[1]]
## [,1] [,2] [,3] [,4]
## [1,] 0 0 0 10.7415
## [2,] 0 0 0 0.0000
## [3,] 0 0 0 0.0000
## [4,] 0 0 0 0.0000
matC(my_comadre)
## [[1]]
## [,1] [,2] [,3] [,4]
## [1,] 0 0 0 0
## [2,] 0 0 0 0
## [3,] 0 0 0 0
## [4,] 0 0 0 0
We can get life history metrics
# the time required for a population to
# increase by a factor of R0 (the net reproductive rate)
gen_time(matU = matU(my_comadre)[[1]], matF(my_comadre)[[1]]) # 17.15
## [1] 17.15455
# the average parent-offspring age difference # 16.22
gen_time(matU = matU(my_comadre)[[1]], matF(my_comadre)[[1]], method = "age_diff")
## [1] 16.22243
# expected age at reproduction for a cohort 5.69
Rage::gen_time(matU = matU(my_comadre)[[1]], matF(my_comadre)[[1]], method = "cohort")
## [1] 5.696541
Rage::life_expect_mean(matU = matU(my_comadre)[[1]], start = 1)
## [1] 1.484286
sum(Rage::life_expect_mean(matU(my_comadre)[[1]], start = NULL))
## [1] 21.87715
Plot the life cycle diagram.
plot_life_cycle(matA = matA(my_comadre)[[1]],
stages = c("eggs/hatchling", "juvenile-early", "juvenile-late", "adult"))
Using survival values from base model.
# The probability of reaching reproductive maturity before death for an egg.
# Expressed as proportion (0 - 1)
# # 0.008000029 # 0.8%
mp_poun <- mature_prob(matU = matU(my_comadre)[[1]], matR = matF(my_comadre)[[1]])
# To produce one reproductive adult need 125 eggs.
1/mp_poun # 125
## [1] 124.9996
# 7 years for a female to lay sufficient eggs
ceiling(125 / 20)
## [1] 7
t0 <- 1000 # number of eggs
n_t1 <- t0 * 0.2 # 200
n_t2 <- n_t1 * 0.333333 # 66.6
n_t3 <- n_t2 * 0.333333 # 22.2
n_t4 <- n_t3 * 0.285714 # 6.3
n_t5 <- n_t4 * 0.285714 # 1.8
n_t6 <- n_t5 * 0.93 # 1.6
n_t7 <- n_t6 * 0.93 # 1.5
n_t8 <- n_t7 * 0.93 # 1.4
n_t9 <- n_t8 * 0.93 # 1.3
ceiling(t0 / n_t5) # 552 eggs to produce one adult (maturity 5 years)
## [1] 552
ceiling(t0 / n_t6) # 593 eggs to produce one adult (maturity 6 years)
## [1] 593
ceiling(t0 / n_t7) # 638 eggs to produce one adult (maturity 7 years)
## [1] 638
ceiling(t0 / n_t8) # 686 eggs to produce one adult (maturity 8 years)
## [1] 686
ceiling(t0 / n_t9) # 737 eggs to produce one adult (maturity 9 years)
## [1] 737
ceiling(552 / 20) # 28 years for a female to lay sufficient eggs (maturity 5 years)
## [1] 28
ceiling(737 / 20) # 37 years for a female to lay sufficient eggs (maturity 9 years)
## [1] 37
# project
dout <- plyr::ddply(dt,
c("species", "type", "first_year","akey"), .fun = pop03_doproj)
dout$arun <- 1
# Model summaries
model_sum <- dout |>
group_by(species, type, first_year, lambda,
gen_time, gen_age_diff, life_exp, life_exp_adult, mat_prob, eggs_to_adult) |>
summarise(fem_t0 = max(fem_t0),
fem_min = min(fem),
fem_max = max(fem)) |>
ungroup()
lambda_n <- length(unique(model_sum$lambda)) # 50
lambda_mean <- mean(model_sum$lambda) # 0.9432
lambda_sd <- sd(model_sum$lambda) # 0.1506
lambda_min <- min(model_sum$lambda) # 0.4659
lambda_max <- max(model_sum$lambda) # 1.1539
# Export for future use
saveRDS(dout, "inst/other/dout.rds")
# Stochastic
#data frame with runs for processing
#nruns <- 100 # 100 gives same pattern as 50
nruns <- 50
dt_stoch <- dt[rep(seq_len(nrow(dt)), nruns), ]
dt_stoch$arun <- rep(1:nruns, each = nrow(dt))
# Approx 90 minutes. 1,212,000 rows. Projections quick. Summaries slow.
# 11:32 - 13:14
dout_stoch <- plyr::ddply(dt_stoch,
c("arun", "species", "type", "first_year","akey"),
.fun = pop03_doproj_stoch)
table(dout_stoch$model)
table(dout_stoch$type)
model_sum_stoch <- dout_stoch |>
group_by(species, type, model, first_year, lambda, lambda_q75,
gen_time, gen_age_diff_med, gen_age_q75,
life_exp_med, life_exp_adult_med, mat_prob_med, mat_prob_q75,
eggs_to_adult_med, eggs_to_adult_q75) |>
summarise(acount = n(),
fem_t0 = max(fem_t0),
fem_min = min(fem),
fem_max = max(fem)) |>
ungroup()
# Export for future use
saveRDS(dout_stoch, "inst/other/dout_stoch.rds")
dout <- readRDS("inst/other/dout.rds")
dout_stoch <- readRDS("inst/other/dout_stoch.rds")
# Combine data for plotting
dout_all <- dplyr::bind_rows(dout |> dplyr::select(arun, model, type, first_year,
akey, ayear,
lambda, gen_time, gen_age_diff,
life_exp, life_exp_adult,
mat_prob, eggs_to_adult,
fem, fem_t0, fem_diff, change50_flag,
change30_flag,
double_flag) |>
dplyr::mutate(lambda_lcl = NA, lambda_ucl = NA,
lambda_sd = NA, gen_sd = NA),
dout_stoch |> dplyr::select(arun, model, type, first_year,
akey, ayear,
lambda, lambda_lcl, lambda_ucl,
lambda_sd, gen_time, gen_sd,
gen_age_diff,
life_exp, life_exp_adult,
mat_prob, eggs_to_adult,
fem, fem_t0, fem_diff, change50_flag,
change30_flag,
double_flag))
# Limit adult female number to maximum (20% above original for baseline).
summary(dout_all$fem)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0 0 1 1389 11 16560701
dout_all[which(dout_all$fem > ceiling_threshold), 'fem' ] <- ceiling_threshold
# summary(dout_all$fem)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.000000 0.005763 1.485889 4.522590 11.240127 12.000000
# Factors in right order
dout_all$modelf <- 1
dout_all[which(dout_all$model=="Stochastic uniform") , 'modelf'] <- 2
dout_all[which(dout_all$model=="Stochastic equal") , 'modelf'] <- 3
dout_all[which(dout_all$model=="Stochastic bad x2") , 'modelf'] <- 4
dout_all[which(dout_all$model=="Stochastic bad x4") , 'modelf'] <- 5
dout_all$modelf <- as.factor(dout_all$modelf)
levels(dout_all$modelf) <- c("Deterministic", "Stochastic uniform",
"Stochastic equal", "Stochastic bad x2",
"Stochastic bad x4")
unique(dout_all$modelf)
table(dout_all$modelf)
dout_all$typef <- 1
dout_all[which(dout_all$type=="female-hunt 2.5%") , 'typef'] <- 2
dout_all[which(dout_all$type=="female-hunt 5%") , 'typef'] <- 3
dout_all[which(dout_all$type=="female-hunt 10%") , 'typef'] <- 4
dout_all[which(dout_all$type=="female-hunt 25%") , 'typef'] <- 5
dout_all[which(dout_all$type=="female-hunt 50%") , 'typef'] <- 6
dout_all$typef <- as.factor(dout_all$typef)
levels(dout_all$typef) <- c("base", "female-hunt 2.5%",
"female-hunt 5%",
"female-hunt 10%", "female-hunt 25%",
"female-hunt 50%")
table(dout_all$typef)
# first year survival
dout_all$first_yearf <- as.factor(dout_all$first_year)
fylev <- paste("first-year\nsurvival\n", seq(0, 0.9, by = 0.1), sep = "")
levels(dout_all$first_yearf) <- fylev
# Export for future use. 17/6/2024 - 1,218,060 rows 21 columns.
saveRDS(dout_all, "inst/other/dout_all.rds")
Model summaries to link with scenarios. Summaries across stochastic runs for each model.
RedList guidelines p38
For example, upper and lower quartiles of the projected magnitude of the
future reduction (i.e., reductions with 25% and 75% probability) may be
considered to represent a plausible range of projected
reduction
Green Status Presence
Present when population is at least 1 as measured by lower quartile
population.
Green Status Viability
Viable when not declining: lower quartile of lambda >= 1. A species
is considered viable in a spatial unit if application of the Regional
Red List Guidelines to the population in the spatial unit would result
in a categorization of ‘Least Concern’ OR ‘Near Threatened and Not
Declining’.
Green Status Functional
Functional when population lower quartile is 10 times base
level.
# load data. 2 August 2024 with 1218060 rows and 21 columns
dout_all <- readRDS("inst/other/dout_all.rds")
# 101 years, 5 model types (deterministic with 6060, stochastic with 6060 * 50),
# 6 harvest levels, 10 first year levels
# Make unique model ID. boot mean is same as mean (at least to 6 decimal places)
# 300 projection models.
model_ref <- dout_all |>
group_by(akey, modelf, typef, first_year, arun, lambda, gen_time,
gen_age_diff, life_exp, life_exp_adult, eggs_to_adult) |>
summarise(yc = length(unique(ayear))) |>
ungroup() |>
group_by(akey, modelf, typef, first_year) |>
summarise(count_runs = length(unique(arun)),
count_years = min(yc),
lambda_mean = mean(lambda),
lambda_min = min(lambda),
lambda_max = max(lambda),
lambda_sd = sd(lambda),
lambda_boot_lcl = Hmisc::smean.cl.boot(lambda)["Lower"],
lambda_boot_ucl = Hmisc::smean.cl.boot(lambda)["Upper"],
lambda_q25 = quantile(lambda, probs = 0.25, na.rm = TRUE),
lambda_q75 = quantile(lambda, probs = 0.75, na.rm = TRUE),
gen_mean = mean(gen_time),
gen_med = median(gen_time),
gen_min = min(gen_time),
gen_max = max(gen_time),
gen_sd = sd(gen_time),
gen_boot_lcl = Hmisc::smean.cl.boot(gen_time)["Lower"],
gen_boot_ucl = Hmisc::smean.cl.boot(gen_time)["Upper"],
gen_q25 = quantile(gen_time, probs = 0.25, na.rm = TRUE),
gen_q75 = quantile(gen_time, probs = 0.75, na.rm = TRUE),
gen_age_mean = mean(gen_age_diff),
gen_age_med = median(gen_age_diff),
gen_age_min = min(gen_age_diff),
gen_age_max = max(gen_age_diff),
gen_age_sd = sd(gen_age_diff),
gen_age_boot_lcl = Hmisc::smean.cl.boot(gen_age_diff)["Lower"],
gen_age_boot_ucl = Hmisc::smean.cl.boot(gen_age_diff)["Upper"],
gen_age_q25 = quantile(gen_age_diff, probs = 0.25, na.rm = TRUE),
gen_age_q75 = quantile(gen_age_diff, probs = 0.75, na.rm = TRUE),
life_exp_mean = mean(life_exp),
life_exp_med = median(life_exp),
life_exp_min = min(life_exp),
life_exp_max = max(life_exp),
life_exp_sd = sd(life_exp),
life_exp_boot_lcl = Hmisc::smean.cl.boot(life_exp)["Lower"],
life_exp_boot_ucl = Hmisc::smean.cl.boot(life_exp)["Upper"],
life_exp_q25 = quantile(life_exp, probs = 0.25, na.rm = TRUE),
life_exp_q75 = quantile(life_exp, probs = 0.75, na.rm = TRUE),
life_exp_adult_mean = mean(life_exp_adult),
life_exp_adult_med = median(life_exp_adult),
life_exp_adult_min = min(life_exp_adult),
life_exp_adult_max = max(life_exp_adult),
life_exp_adult_sd = sd(life_exp_adult),
life_exp_adult_boot_lcl = Hmisc::smean.cl.boot(life_exp_adult)["Lower"],
life_exp_adult_boot_ucl = Hmisc::smean.cl.boot(life_exp_adult)["Upper"],
life_exp_adult_q25 = quantile(life_exp_adult, probs = 0.25, na.rm = TRUE),
life_exp_adult_q75 = quantile(life_exp_adult, probs = 0.75, na.rm = TRUE),
eggs_to_adult_mean = mean(eggs_to_adult),
eggs_to_adult_med = median(eggs_to_adult),
eggs_to_adult_min = min(eggs_to_adult),
eggs_to_adult_max = max(eggs_to_adult),
eggs_to_adult_sd = sd(eggs_to_adult),
eggs_to_adult_boot_lcl = Hmisc::smean.cl.boot(eggs_to_adult)["Lower"],
eggs_to_adult_boot_ucl = Hmisc::smean.cl.boot(eggs_to_adult)["Upper"],
eggs_to_adult_q25 = quantile(eggs_to_adult, probs = 0.25, na.rm = TRUE),
eggs_to_adult_q75 = quantile(eggs_to_adult, probs = 0.75, na.rm = TRUE)
) |>
ungroup() |>
arrange(typef, first_year, modelf) |>
mutate(modelid = paste(akey, as.numeric(modelf), sep = "_"),
modelkey = row_number()) |>
relocate(modelkey, modelid)
# Estimates of generation time for summary.
# Bienvenu and Legendre (2015) https://doi.org/10.1086%2F681104 -
# defined as the mean age of mothers at birth.
# IUCN definition - is the mean age at which a cohort of individuals produce offspring.
gen_mean <- model_ref |>
filter(lambda_q75 >= 1, first_year < 0.5) |> pull(gen_age_mean) |> mean()
gen_median <- model_ref |>
filter(lambda_q75 >= 1, first_year < 0.5) |> pull(gen_age_mean) |> median()
gen_q25 <- model_ref |>
filter(lambda_q75 >= 1, first_year < 0.5) |> pull(gen_age_mean) |>
quantile(probs = 0.25)
gen_q75 <- model_ref |>
filter(lambda_q75 >= 1, first_year < 0.5) |> pull(gen_age_mean) |>
quantile(probs = 0.75)
gen_3 <- gen_mean * 3
gen_3_q25 <- gen_q25 * 3
gen_3_q75 <- gen_q75 * 3
# Add population changes.
model_ref2 <- model_ref |> left_join(dout_all |>
filter(ayear %in% c(ceiling(gen_3), ceiling(gen_3_q25), ceiling(gen_3_q75), 100)) |>
group_by(akey, modelf, typef, first_year, ayear, arun) |>
summarise(fem_t0 = mean(fem_t0),
femnew = mean(fem),
femnew_diff = mean(fem_diff)
) |>
ungroup() |>
group_by(akey, modelf, typef, first_year, ayear) |>
summarise(fem_t0 = mean(fem_t0),
fem = mean(femnew),
fem_sd = sd(femnew, na.rm = TRUE),
fem_q25 = quantile(femnew, probs = 0.25, na.rm = TRUE),
fem_q75 = quantile(femnew, probs = 0.75, na.rm = TRUE),
fem_diff = mean(femnew_diff),
fem_diff_sd = sd(femnew_diff, na.rm = TRUE),
fem_diff_q25 = quantile(femnew_diff, probs = 0.25, na.rm = TRUE),
fem_diff_q75 = quantile(femnew_diff, probs = 0.75, na.rm = TRUE)) |>
ungroup() |>
mutate(ayear = paste("t", ayear, sep="")) |>
pivot_wider(names_from = ayear,
values_from = c(fem, fem_sd, fem_q25, fem_q75, fem_diff,
fem_diff_sd, fem_diff_q25, fem_diff_q75)) |>
arrange(typef, first_year, modelf)
)
# Add time to change by 50% and 30%
model_ref3 <- model_ref2 |>
left_join(dout_all |>
filter(change50_flag == 1) |>
group_by(akey, modelf, typef, first_year, arun,
change50_flag) |>
summarise(acount = n(),
change50_yf = min(ayear),
change50_yl = max(ayear)) |>
ungroup() |>
group_by(akey, modelf, typef, first_year) |>
summarise(change50_ymean = floor(mean(change50_yf)),
change50_ysd = sd(change50_yf),
change50_yq25 = floor(quantile(change50_yf, probs = 0.25, na.rm = TRUE)),
change50_yq75 = floor(quantile(change50_yf, probs = 0.75, na.rm = TRUE))) |>
ungroup()
) |> left_join(dout_all |>
filter(change30_flag == 1) |>
group_by(akey, modelf, typef, first_year, arun,
change30_flag) |>
summarise(acount = n(),
change30_yf = min(ayear),
change30_yl = max(ayear)) |>
ungroup() |>
group_by(akey, modelf, typef, first_year) |>
summarise(change30_ymean = floor(mean(change30_yf)),
change30_ysd = sd(change30_yf),
change30_yq25 = floor(quantile(change30_yf, probs = 0.25, na.rm = TRUE)),
change30_yq75 = floor(quantile(change30_yf, probs = 0.75, na.rm = TRUE))) |>
ungroup()
)
# Export for future use
saveRDS(model_ref3, "inst/other/model_lookup.rds")
write.csv2(model_ref3, "inst/other/model_lookup.csv", row.names = FALSE)
# summaries of 10050 runs.
#model_all_sum <- dout_all |>
# group_by(akey, modelf, typef, first_year, arun) |>
# summarise(count_years = length(ayear),
# count_years_unique = length(unique(ayear))) |>
# ungroup()
Compare lambda. Use this as the basis to establish number of runs to use.
model_lookup <- readRDS("inst/other/model_lookup.rds")
model_lookup |>
mutate(recover_flag = if_else(fem_diff_t100 > 0, 1, 0)) |>
group_by(recover_flag) |>
summarise(model_count = n(),
lambda_avg = mean(lambda_mean),
lambda_min = min(lambda_mean),
lambda_max = max(lambda_mean))
## # A tibble: 2 × 5
## recover_flag model_count lambda_avg lambda_min lambda_max
## <dbl> <int> <dbl> <dbl> <dbl>
## 1 0 212 0.842 0.466 0.995
## 2 1 88 1.04 0.973 1.15
# 18 where mean lambda is less than one but population increases.
# 100, 144 and 198 with lowest lambda. Results uncertain
model_lookup |>
filter(fem_diff_t100 > 0, lambda_mean < 1) |>
arrange(lambda_mean)
## # A tibble: 19 × 102
## modelkey modelid akey modelf typef first_year count_runs count_years
## <int> <chr> <chr> <fct> <fct> <dbl> <int> <int>
## 1 193 poun_female-hu… poun… Stoch… fema… 0.8 50 101
## 2 144 poun_female-hu… poun… Stoch… fema… 0.8 50 101
## 3 100 poun_female-hu… poun… Stoch… fema… 0.9 50 101
## 4 198 poun_female-hu… poun… Stoch… fema… 0.9 50 101
## 5 149 poun_female-hu… poun… Stoch… fema… 0.9 50 101
## 6 35 poun_base_0.6_5 poun… Stoch… base 0.6 50 101
## 7 84 poun_female-hu… poun… Stoch… fema… 0.6 50 101
## 8 73 poun_female-hu… poun… Stoch… fema… 0.4 50 101
## 9 40 poun_base_0.7_5 poun… Stoch… base 0.7 50 101
## 10 133 poun_female-hu… poun… Stoch… fema… 0.6 50 101
## 11 24 poun_base_0.4_4 poun… Stoch… base 0.4 50 101
## 12 89 poun_female-hu… poun… Stoch… fema… 0.7 50 101
## 13 45 poun_base_0.8_5 poun… Stoch… base 0.8 50 101
## 14 50 poun_base_0.9_5 poun… Stoch… base 0.9 50 101
## 15 94 poun_female-hu… poun… Stoch… fema… 0.8 50 101
## 16 18 poun_base_0.3_3 poun… Stoch… base 0.3 50 101
## 17 29 poun_base_0.5_4 poun… Stoch… base 0.5 50 101
## 18 78 poun_female-hu… poun… Stoch… fema… 0.5 50 101
## 19 138 poun_female-hu… poun… Stoch… fema… 0.7 50 101
## # ℹ 94 more variables: lambda_mean <dbl>, lambda_min <dbl>, lambda_max <dbl>,
## # lambda_sd <dbl>, lambda_boot_lcl <dbl>, lambda_boot_ucl <dbl>,
## # lambda_q25 <dbl>, lambda_q75 <dbl>, gen_mean <dbl>, gen_med <dbl>,
## # gen_min <dbl>, gen_max <dbl>, gen_sd <dbl>, gen_boot_lcl <dbl>,
## # gen_boot_ucl <dbl>, gen_q25 <dbl>, gen_q75 <dbl>, gen_age_mean <dbl>,
## # gen_age_med <dbl>, gen_age_min <dbl>, gen_age_max <dbl>, gen_age_sd <dbl>,
## # gen_age_boot_lcl <dbl>, gen_age_boot_ucl <dbl>, gen_age_q25 <dbl>, …
# 10 where t100 females can be more or less than t0.
model_lookup |>
mutate(flag_uncertain = if_else((fem_q25_t100 < fem_t0) & (fem_q75_t100 > fem_t0),
1, 0)) |>
filter(flag_uncertain == 1) |>
arrange(lambda_mean)
## # A tibble: 9 × 103
## modelkey modelid akey modelf typef first_year count_runs count_years
## <int> <chr> <chr> <fct> <fct> <dbl> <int> <int>
## 1 193 poun_female-hun… poun… Stoch… fema… 0.8 50 101
## 2 144 poun_female-hun… poun… Stoch… fema… 0.8 50 101
## 3 100 poun_female-hun… poun… Stoch… fema… 0.9 50 101
## 4 128 poun_female-hun… poun… Stoch… fema… 0.5 50 101
## 5 35 poun_base_0.6_5 poun… Stoch… base 0.6 50 101
## 6 84 poun_female-hun… poun… Stoch… fema… 0.6 50 101
## 7 73 poun_female-hun… poun… Stoch… fema… 0.4 50 101
## 8 40 poun_base_0.7_5 poun… Stoch… base 0.7 50 101
## 9 24 poun_base_0.4_4 poun… Stoch… base 0.4 50 101
## # ℹ 95 more variables: lambda_mean <dbl>, lambda_min <dbl>, lambda_max <dbl>,
## # lambda_sd <dbl>, lambda_boot_lcl <dbl>, lambda_boot_ucl <dbl>,
## # lambda_q25 <dbl>, lambda_q75 <dbl>, gen_mean <dbl>, gen_med <dbl>,
## # gen_min <dbl>, gen_max <dbl>, gen_sd <dbl>, gen_boot_lcl <dbl>,
## # gen_boot_ucl <dbl>, gen_q25 <dbl>, gen_q75 <dbl>, gen_age_mean <dbl>,
## # gen_age_med <dbl>, gen_age_min <dbl>, gen_age_max <dbl>, gen_age_sd <dbl>,
## # gen_age_boot_lcl <dbl>, gen_age_boot_ucl <dbl>, gen_age_q25 <dbl>, …
mean_values <- model_ref |>
group_by(typef) |>
summarise(l_mean = mean(lambda_mean),
l_lcl = Hmisc::smean.cl.boot(lambda_mean)["Lower"],
l_ucl = Hmisc::smean.cl.boot(lambda_mean)["Upper"],
)
fig_model_lambda <- model_ref |>
ggplot(aes(x = first_year, y = lambda_mean, colour = modelf)) +
geom_hline(data = mean_values, aes(yintercept = l_mean),
colour = "black", linewidth = 0.6,
linetype = "dashed") +
geom_hline(data = mean_values, aes(yintercept = l_lcl),
colour = "grey", linewidth = 0.6,
linetype = "dashed") +
geom_hline(data = mean_values, aes(yintercept = l_ucl),
colour = "grey", linewidth = 0.6,
linetype = "dashed") +
geom_hline(yintercept = 1,
colour = "magenta", linewidth = 0.4) +
geom_point() +
stat_smooth(se = FALSE) +
scale_x_continuous("first year survival and graduation",
breaks = c(0, 0.2, 0.4, 0.6, 0.8)) +
scale_color_viridis_d() +
facet_wrap(~ typef, ncol = 6) +
theme(legend.position="top") +
guides(colour=guide_legend(nrow=2,byrow=TRUE)) +
labs(y = "lambda", colour = "model")
# save
png(file = "inst/other/fig_unifilis_lambda.png", bg = "transparent",
type = c("cairo"),
width = 9, height = 4, units = "in", res=600)
fig_model_lambda
invisible(dev.off())
Check plot.
Model lambda for Podocnemis unifilis
fig_proj <- dout_all |>
ggplot(aes(x = ayear, y = fem, colour = modelf)) +
geom_point(size=0.1, alpha=0.2) +
stat_smooth(se = FALSE) +
scale_colour_viridis_d("model") +
scale_y_continuous(limits = c(0, ceiling_threshold),
breaks = seq(2, 12, by = 2)) +
labs(y="Reproductive females",
x="Time (years)",
) +
theme_bw() +
theme(plot.title.position = "plot") +
facet_grid(first_yearf ~ typef)
# save
png(file = "inst/other/fig_unifilis_projections.png", bg = "transparent",
type = c("cairo"),
width = 10, height = 8, units = "in", res=600)
fig_proj
invisible(dev.off())
Check plot.